Non-Bayesian particle filters 



Alexandre J. Chorin and Xuemin Tu 

Department of Mathematics, University of California at Berkeley 

and 

Lawrence Berkeley National Laboratory 
Berkeley, CA, 94720 

Abstract 

Particle filters for data assimilation in nonlinear problems use "par- 
ticles" (replicas of the underlying system) to generate a sequence of 
probability density functions (pdfs) through a Bayesian process. This 
can be expensive because a significant number of particles has to be 
used to maintain accuracy. We offer here an alternative, in which 
the relevant pdfs are sampled directly by an iteration. An example is 
discussed in detail. 
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1 Introduction. 

There are many problems in science in which the state of a system must 
be identified from an uncertain equation supplemented by a stream of noisy 
data (see e.g. pQ). A natural model of this situation consists of a stochastic 
differential equation (SDE): 

<ix = f (x, t) dt + g(x, t) dw, (1) 

where x = (x±, X2, ■ ■ ■ , x m ) is an m-dimensional vector, dw is m-dimensional 
Brownian motion, f is an m-dimensional vector function, and g is a scalar 
(i.e., an m by m diagonal matrix of the form gl, where g is a scalar and / is 
the identity matrix). The Brownian motion encapsulates all the uncertainty 
in this equation. The initial state x(0) is assumed given and may be random 
as well. 

As the experiment unfolds, it is observed, and the values b n of a mea- 
surement process are recorded at times t n ; for simplicity assume t n = n5, 
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where 6 is a fixed time interval and n is an integer. The measurements are 
related to the evolving state x(t) by 



where h is a fc-dimensional, generally nonlinear, vector function with k < m, 
G is a diagonal matrix, x n = x(n<5), and W n is a vector whose components 
are independent Gaussian variables of mean and variance 1, independent 
also of the Brownian motion in equation ([1]). The task is to estimate x on 
the basis of equation ([T]) and the observations (J2]). 

If the system ([T]) is linear and the data are Gaussian, the solution can be 
found via the Kalman-Bucy filter. In the general case, it is natural to try to 
estimate x as the mean of its evolving probability density. The initial state 
x is known and so is its probability density; all one has to do is evaluate 
sequentially the density P n+ i of x n+1 given the probability density P n of x n 
and the data b n+1 . This can be done by following "particles" (replicas of 
the system) whose empirical distribution approximates P n . In a Bayesian 
filter (see e.g El H El El [3 El [9] , one uses the pdf P n and equation $TJ to 
generate a prior density, and then one uses the new data b n+1 to generate 
a posterior density P n +\- In addition, one may have to sample backward to 
take into account the information each measurement provides about the past 
and avoid having too many identical particles. Evolving particles is typically 
expensive, and the backward sampling, usually done by Markov chain Monte 
Carlo (MCMC), can be expensive as well, because the number of particles 
needed can grow catastrophically (see e.g. [TO]). 

In this paper we offer an alternative to the standard approach, in which 
P n +i is sampled directly without recourse to Bayes' theorem and backward 
sampling, if needed, is done by chainless Monte Carlo [11] . Our direct sam- 
pling is based on a representation of a variable with density P n+ i by a col- 
lection of functions of Gaussian variables parametrized by the support of P n , 
with parameters found by iteration. The construction is related to chainless 
sampling as described in [TTJ. The idea in chainless sampling is to produce 
a sample of a large set of variables by sequentially sampling a growing se- 
quence of nested conditionally independent subsets. As observed in [121 H3] , 
chainless sampling for a SDE reduces to interpolatory sampling, as explained 
below. Our construction will be explained in the following sections through 
an example where the position of a ship is deduced from the measurements 
of an azimuth, already used as a test bed in [61 [HI [15] . 
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2 Sampling by interpolation and iteration. 



First we explain how to sample via interpolation and iteration in a simple 
example, related to the example and the construction in [12]. Consider the 
scalar SDE 

dx = f(x,t)dt + y/a dw; (3) 

we want to find sample paths x = x(t), < t < 1, subject to the conditions 
x(0) = 0,x(l) = X. 

Let N(a,v) denote a Gaussian variable with mean a and variance v. We 
first discretize equation ([31) on a regular mesh t ,* 1 , . . . ,t N , where t n = n5, 
5 = 1/N, < n < N, with x n = x(t n ), and, following [12], use a balanced 
implicit discretization [T6l [T7]: 

x n+l =x n + y^n f n^ + _ x ")f( x ")S + W n+ \ 

where f'(x n ,t n ) = J^(x n ,t n ) and W n+1 is N(0,cr/N). The joint probability 
density of the variables x 1 , . . . , x^ -1 is Z~ x exp(— Y2o ^™)> where Z is the 
normalization constant and 

^ _ ((l-5f')(x n+l -x n )-5f) 2 
2a5 

(x n+1 -x n -Sf/{l -8f')f 
= 2a~ n ' 

where /, /' are functions of the x\ and a n = a8/(l — Sf') 2 (see [IE])- One 
can obtain sample solutions by sampling this density, e.g. by MCMC, or one 
can obtain them by interpolation (chainless sampling), as follows. 

Consider first the special case f(x,t) = f(t), so that in particular /' = 0. 
Each increment x n+l — x n is now a N(a n , cr/N) variable, with the a n = f(t n )S 
known explicitly. Let iV be a power of 2. Consider the variable x N ^ 2 . On 
one hand, 

N/2 

x N/2 = ^(x- - a;™" 1 ) = N(Ai, V x ), 

1 

where A\ = Yli^ 2 a m V± = a/2. On the other hand, 

N 

x = x N / 2 + Yl o^-^ -1 ), 

7V/2+1 



3 



so that 

z*' 2 = N(A 2 ,V 2 ), 

with 

N-l 

A 2 = X - a n, V2 = Vl 

N/2+1 

The pdf of x N l 2 is the product of the two pdfs; one can check that 

(x-A 1 ) 2 \ / (x-A 2 ) 2 



exp I — — exp 



2V X J 1 \ 2V 2 
exp ( — — I exp(-0), 

where v = jgfe, a = Vl ^+^ 2 , and = e"* is the probability of 

getting from the origin to X, up to a normalization constant. 

Pick a sample £1 from the N(0, 1) density; one obtains a sample of x N l 2 
by setting i N / 2 = a + y^i- Given a sample of i^ 2 one can similarly sample 
i w / 4 , x 3Ar / 4 , then s^/ 8 , x 3Ar//s , etc., until all the x^ have been sampled. If 
we define £ = (£1,62, • • • , £n-i), then for each choice of £ we find a sample 
(x 1 , . . . , i^" 1 ) such that 

-.1 1 1 C2 



exp ( — — - - ) oxp 

exp 



2 / \ 2d 

(rr 1 — x° — a ) 2 (x 2 — x 1 — ax) 2 



2a /N 2a /N 

(x N -x N ~ l -a N ^) 2 ' 
2a /N 



(4) 



where the factor exp ^— — j on the left is the probability of the fixed 

end value X up to a normalization constant. In this linear problem, this 
factor is the same for all the samples and therefore harmless. One can repeat 
this sampling process for multiple choices of the variables each sample of 
the corresponding set of x n is independent of any previous samples of this 
set. 

Now return to the general case. The functions /, /' are now functions of 
the xK We obtain a sample of the probability density we want by iteration. 
First pick S = (£1,^2, • • • , 6v-i), where each £j is drawn independently from 
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the iV(0, 1) density (this vector remains fixed during the iteration). Make 
a first guess x° = (xJ,Xq, . . . ,Xq _1 ) (for example, if X ^ 0, pick x = 0). 
Evaluate the functions /, /' at x J (note that now /' ^ 0, and therefore the 
variances of the various increments are no longer constants). We are back in 
previous case, and can find values of the increments — x™ +1 corresponding 
to the values of /, /' we have. Repeat the process starting with the new 
iterate. If the vectors x J converge to a vector x = (x 1 , . . . , x^ -1 ), we obtain, 
in the limit, equation (j3J), where now on the right side a depends on n so 
that a = a n , and both a n , a n are functions of the final x. The left hand side 
of P| becomes: 

exp exp 1 



2 y 



Note that now the factor exp ^— ^^—^ ) is different from sample to sam- 
ple, and changes the relative weights of the different samples. In averaging, 
one should take this factor as weight, or resample as described at the end of 
the following section. In order to obtain more uniform weights, one also can 
use the strategies in [TTJ [12] . 

One can readily see that the iteration converges if KTM < 1, where K 
is the Lipshitz constant of /, T is the length of the interval on which one 
works (here T = 1), and M is the maximum norm of the vectors x J+1 — 
x 3 '. If this inequality is not satisfied for the iteration above, it can be re- 
established by a suitable underrelaxation. One should course choose N large 
enough so that the results are converged in N. We do not provide more 
details here because they are extraneous to our purpose, which is to explain 
chainless / interpolatory sampling and the use of reference variables in a simple 
context. 



3 The ship azimuth problem. 

The problem we focus on is discussed in [HI HH US], where it is used to 
demonstrate the capabilities of particular Bayesian filters. A ship sets out 
from a point (xo,yo) in the plane and undergoes a random walk, 

x n+l =x n + dx n+l^ 

y n+1 =y n + dy n+1 , (5) 
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for n > 0, and with x° = y° given, and dx n+l = N(dx n ,a), dy n+1 = 
N(dy n , a), i.e., each displacement is a sample of a Gaussian random variable 
whose variance a does not change from step to step and whose mean is the 
value of the previous displacement. An observer makes noisy measurements 
of the azimuth arctan(y ri /x n ), recording 

b n = arctan^ + N(0,s). (6) 

where the variance s is also fixed; here the observed quantity b is scalar and 
is not be denoted by a boldfaced letter. The problem is to reconstruct the 
positions x n = (x n , y n ) from equations (151161) . We take the same parameters as 
[6]: x = 0.01, y = 20, dx 1 = 0.002, dy 1 = -0.06, a = 1 ■ lO" 6 , s = 25 ■ 10" 6 . 
We follow numerically M particles, all starting from X 4 ° = xq, Y® = yo, 
as described in the following sections, and we estimate the ship's position 
at time n5 as the mean of the locations X" = (X™, Y™), i — 1, . . . , M of the 
particles at that time. The authors of [6] also show numerical results for runs 
with varying data and constants; we discuss those refinements in section 6 
below. 



4 Forward step. 

Assume we have a collection of M particles X n at time t n = nS whose 
empirical density approximates P n ; now we find increments c£X n+1 such that 
the empirical density of X" +1 = X n + dX. n+1 approximates P n+ i- P n +i is 
known implicitly: it is the product of the density that can be deduced from 
the SDE and the one that comes from the observations, with the appropriate 
normalization. If the increments were known, their probability p (the density 
P n +\ evaluated at the resulting positions X n+1 ) would be known, so p is a 
function of dX n+1 , p = p(dX. n+1 ). For each particle i, we are going to sample 
a Gaussian reference density, obtain a sample of probability p, then solve (by 
iteration) the equation 

p = pidxr 1 ) (?) 

to obtain dX™ +1 . 

Define f(x,y) = arctan(|//x) and f n = f(X n ,Y n ). We are working on 
one particle at a time, so the index i can be temporarily suppressed. Pick 
two independent samples £ y from a N(0, 1) density (the reference density 
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in the present calculation), and set p = ^ exp (^—^f — yj; the variables 
£ x , £ y remain unchanged until the end of the iteration. We are looking for 
displacements dX n+1 , dY n+l , and parameters a x , a y ,v x ,v y , (ft, such that: 

/ (dX n+1 - dX n ) 2 (dY n+1 -dY n ) 2 
2vrp=exp [ - — 



2a 2a 

(/ n+1 — } ) n + 1 \ 2 ^ 



2s 



exp(0) 



^expf-^ +1 -^ 2 -^ +1 -^H (8) 
V 2v x 2v y J y 1 

The first equality states what we wish to accomplish: find increments dX n+1 , 
dY n+1 , functions respectively of £, x ,£, y , whose probability with respect to 
P n+ i is p. The factor is needed to normalize this term (0 is called below a 
"phase"). The second equality says how the goal is reached: we are looking 
for parameters a x ,a y ,v x ,v y , (all functions of X n ) such that the increments 
are samples of Gaussian variables with these parameters, with the assumed 
probability. One should remember that in our example the mean of dX n+1 
is dX n , and similarly for dY n+1 . We are not representing P n+ \ as a function 
of a single Gaussian- there is a different Gaussian for every value of X n . 

To satisfy the second equality we set up an iteration for vectors dX n+1 ' J '(= 
g?X j for brevity) that converges to <£X n+1 . Start with <£X° = 0. We now 
explain how to compute dX.i +1 given cOC . 

Approximate the observation equation (jB]) by 

/(X') + f x ■ (dXi +1 - dX>) + /„ • (dYi +1 - dYi) = b n+1 + N(0, s), (9) 

where the derivatives f x , f y are, like /, evaluated at X J = X n + dXP, i.e., ap- 
proximate the observation equation by its Taylor series expansion around the 
previous iterate. Define a variable rf +x = (f x ■ dX^ +1 + f y ■ dY j+1 ) / f x + f y - 
The approximate observation equation says that r]^ +1 is a N(ai,Vi) variable, 
with 

f-f x - dXi - f y ■ dY? - b n+1 
a\ = 



J X ' J V 



'V 
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On the other hand, from the equations of motion one finds that r]^ +l is 
N(a 2 , v 2 ), with a 2 = (f x ■ dX n + f y ■ dY n ) / ^/ ft + and v 2 = a. Hence the 
pdf of rj^ 1 is, up to normalization factors, 

(x — ai) 2 (x — a 2 ) 2 \ f (x — a^ 2 



" X1>1 2^ 2^J =eXP { 2^» <,X1)! 

where v = a = 6 = t^TK = 

We can also define a variable 7/^ +1 that is a linear combination of dX^ +1 , 
dYi +l and is uncorrelated with if +1 : 

i+ i -A • + /* • 
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The observations do not affect 77+" , so its mean and variance are known. 
Given the means and variances of r]+ +1 one can easily invert the orthog- 
onal matrix that connects them to dX^ +1 , dY^ +1 and find the means and 
variances a x ,v x of dX^ +1 and a y ,v y of dY^ +1 after their modification by the 
observation (the subscripts on a, v are labels, not differentiations). Now one 
can produce values for dX^ +1 , dYi +l : 

dX J+1 = a x + y/%^ x , dY J+1 = a y + y/Vy~£ y , 

where £ x , £ y are the samples from N(0, 1) chosen at the beginning of the 
iteration. This completes the iteration. 

This iteration converges to X n+1 such that /(X n+1 ) = b n+1 + N(0,s), 
and the phases (fP converge to a limit <fi — 0j, where the particle index i 
has been restored. The time interval over which the solution is updated 
in each step is short, and we do not expect any problem with convergence, 
either here or in the next section, and indeed there is none; in all cases the 
iteration converges in a small number of steps. Note that after the iteration 
the variables X" +1 ,y. n+1 are no longer independent- the observation creates 
a relation between them. 

Do this for all the particles. The particles are now samples of P n +i, but 
they have been obtained by sampling different densities (remember that the 
parameters in the Gaussians in equation (jSJ) vary) . One can get rid of this het- 
erogeneity by viewing the factors exp(— <p) as weights and resampling, i.e., for 
each of M random numbers 9k, k — 1, . . . , M drawn from the uniform distri- 
bution on [0, 1], choose a new X£ +1 = X™ +1 such that Z^ 1 Y^j=i ex P( — 4>j) < 
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@k < Z _1 Y^]=i ex P(~ 4>j) (where Z = Y^jLi exp(— (f>j)), and then suppress 
the hat. We have traded the resampling of Bayesian filters for a resampling 
based on the normalizing factors of the several Gaussian densities; this is a 
worthwhile trade because in a Bayesian filter one gets a set of samples many 
of which may have low probability with respect to P n+ i, and here we have a 
set of samples each one of which has high probability with respect to a pdf 
close to P n+ i- 

Note also that the resampling does not have to be done at every step- for 
example, one can add up the phases for a given particle and resample only 
when the ratio of the largest cumulative weight exp(— ^ <&) to the smallest 
such weight exceeds some limit L (the summation is over the weights accrued 
to a particular particle % since the last resampling). If one is worried by 
too many particles being close to each other ("depletion" in the Bayesian 
terminology), one can divide the set of particles into subsets of small size 
and resample only inside those subsets, creating a greater diversity. As will 
be seen in section 6, none of these strategies will be used here and we will 
resample fully at every step. 

5 Backward sampling. 

The algorithm of the previous section is sufficient to create a filter, but 
accuracy may require an additional refinement. Every observation provides 
information not only about the future but also about the past- it may, for 
example, tag as improbable earlier states that had seemed probable before 
the observation was made; one may have to go back and correct the past after 
every observation (this backward sampling is often misleadingly motivated 
solely by the need to create greater diversity among the particles in a Bayesian 
filter). As will be seen below, this backward sampling does not provide a 
significant boost to accuracy in the present problem, but it is described here 
for the sake a completeness. 

Given a set of particles at time (n + 1)5, after a forward step and maybe 
a subsequent resampling, one can figure out where each particle i was in the 
previous two steps, and have a partial history for each particle v. X™" 1 , X™, X™ 
(if resamples had occurred, some parts of that history may be shared among 
several current particles). Knowing the first and the last member of this 
sequence, one can interpolate for the middle term as in section 2, thus pro- 
jecting information backward. This requires that one recompute <£K ra . 
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Let dX. tot = dX. n + dX. n+ ; in the present section this quantity is assumed 
known and remains fixed. In the azimuth problem discussed here, one has 
to deal with the slight complication due to the fact that the mean of each 
increment is the value of the previous one, so that two successive increments 
are related in a slightly more complicated way than usual. The displacement 
dX n is a N(dX n ^ 1 , a) variable, and dX n+l is a N(dX n ,a) variable, so that 
one goes from X n ~ l to X n+1 by sampling first a (2dX n ~ 1 , 4er) variable that 
takes us from X n ~ x to an intermediate point P, with a correction by the 
observation half way up this first leg, and then one samples a N(dX tot , a) 
variable to reach X n+l , and similarly for Y. Let the variable that connects 



X"" 1 to P be dX ncw , so that what replaces dX n is dX new /2. Accordingly, 
we are looking for a new displacement <iX new = (c?X new , dY new ), and for 
parameters < ew , < ew , < ew , < ew such that 



where f n = /(I"" 1 + dX ncw /2, Y n ~ l + dY ncw /2) and are independent 

N(0, 1) Gaussian variables. As in equation (jSj), the first equality embodies 
what we wish to accomplish- find increments, functions of the reference vari- 
ables, that sample the new pdf at time nS defined by the forward motion, 
the constraint imposed by the observation, and by knowledge of the position 
at time (n + l)St. The second equality states that this is done by finding 
particle-dependent parameters for a Gaussian density. 

We again find these parameters as well as the increments by iteration. 
Much of the work is separate for the X and Y components of the equations of 
motion, so we write some of the equations for the X component only. Again 
set up an iteration for variables dJ newJ = dX^ which converge to dX new . 
Start with dX° = 0. To find dX^ +l given dX\ approximate the observation 
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equation flS]), as before, by equation flU]); define again variables rf +1 , r/+ +1 , one 
in the direction of the approximate constraint and one orthogonal to it; in 
the direction of the constraint multiply the pdfs as in the previous section; 
construct new means a^.a^ and new variances v^Vy for dX,dY at time n, 
taking into account the observation at time n, again as before. This also 
produces a phase = 4> . 

Now take into account that the location of the boat at time n+1 is known; 
this creates a new mean a x , a new variance v x , and a new phase <p x , by v = 
^SS 3 , 0* = where a, = 2a\v x = 4t£, a 2 = X tot , v 2 = 



1)1+1)2 ' D1+D2 ' ' X 1)1+1)2 

a. Finally, find a new interpolated position dX^ +1 = a" cw /2 + y/v™ w £ x (the 
calculation for dY^ +l is similar, with a phase fa), and we are done. The total 
phase for in this iteration is (j) = fa + (f) x + fa. As the iterates dX^ converge 
to dX new , the phases converge to a limit = fa. The probability of a particle 
arriving at the given position at time {n + l)St having been determined in 
the forward step, there is no need to resample before comparing samples. 
Once one has the values of X new , a forward step gives corrected values of 
X n+1 ; one can use this interpolation process to correct estimates of X fc by 
subsequent observations for k — n — 1, k — n — 2, . . . , as many as are useful. 



6 Numerical results. 

Before presenting examples of numerical results for the azimuth problem, we 
discuss the accuracy one can expect. A single set of observations for our 
problem relies on 160 samples of a N(0, a) variable. The maximum likeli- 
hood estimate of a given these samples is a random variable with mean a 
and standard deviation .Wo. We estimate the uncertainty in the position 
of the boat by picking a set of observations, then making multiple runs of 
the boat where the random components of the motion in the direction of the 
constraint are frozen while the ones orthogonal to it are sampled over and 
over from the suitable Gaussian density, then computing the distances to the 
fixed observations, estimating the standard deviation of these differences, 
and accepting the trajectory if the estimated standard deviation is within 
one standard deviation of the nominal value of s. This process generates a 
family of boat trajectories compatible with the given observations. In Table 
I we display the standard deviations of the differences between the resulting 
paths and the original path that produced the observations after the num- 
ber of steps indicated there (the means of these differences are statistically 
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indistinguishable from zero). This Table provides an estimate of the accu- 
racy we can expect. It is fair to assume that these standard deviations are 
underestimates of the uncertainty- a variation of a single standard deviation 
in s is a strict constraint, and we allowed no variation in a. 

Table I 

Intrinsic uncertainty in the azimuth problem 



step 


x component 


y component 


40 


.0005 


.21 


80 


.004 


.58 


120 


.010 


.88 


160 


.017 


.95 



If one wants reliable information about the performance of the filter, it 
is not sufficient to run the boat once, record observations, and then use the 
filter to reconstruct the boat's path, because the difference between the true 
path and the reconstruction is a random variable which may be accidentally 
atypically small or atypically large. We have therefore run a large number 
of such reconstructions and computed the means and standard deviations 
of the discrepancies between path and reconstruction as a function of the 
number of steps and of other parameters. In Tables II and III we display the 
means and standard deviations of these discrepancies (not of their mean!) 
in the the x and y components of the paths with 2000 runs, at the steps 
and numbers of particles indicated, with no backward sampling. (Ref. [6] 
used 100 particles). On the average the error is zero, and the error that can 
be expected in any one run is of the order of magnitude of the unavoidable 
error. The standard deviation of the discrepancy is not significantly smaller 
with 2 particles that with 100- the main source of the discrepancy is the 
uncertainty in the data. Most of time one single particle (no resampling) is 
enough; however, a single particle may temporarily stray into low-probability 
areas and creates large arguments and numerical difficulties in the various 
functions used in the program. Two particles with resampling keep each 
other within bounds, because if one of them strays it gets replaced by a 
replica of the other. The various more sophisticated resampling strategies 
at the end of section 4 make no discernible difference here, and backward 
sampling does not help much either, because they too are unable to remedy 
the limitation of the data set. 
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Table Ha 

Mean and standard variation of the discrepancy between synthetic data and 
their reconstruction, 2000 runs, no back step, 100 particles 



n. of steps 


x component 


y component 




mean 


s.d. 


mean 


s.d. 


40 


.0004 


.04 


.0001 


.17 


80 


-.001 


.04 


-.01 


.54 


120 


-.0008 


.07 


-.03 


1.02 


160 


-.002 


.18 


-.05 


1.56 



Table lib 

Mean and standard variation of the discrepancy between synthetic data and 
their reconstruction, 2000 runs, no back step, 2 particles 



n. of steps 


x com 


3onent 


y component 




mean 


s.d. 


mean 


s.d. 


40 


.002 


.17 


-.0004 


.20 


80 


.01 


.43 


-.0006 


.58 


120 


.01 


.57 


.009 


1.08 


160 


.006 


.54 


.01 


1.67 



In Figure 1 we plot a sample boat path, its reconstruction, and the re- 
constructions obtained (i) when the initial data for the reconstruction are 
strongly perturbed (here, the initial data for x, y were perturbed initially by, 
respectively, .1 and .4), and (ii) when the value of a assumed in the recon- 
struction is random: a = N(ao, ecr ), where <7o is the constant value used until 
now and e = 0.4 but the calculation is otherwise identical. This produces 
variations in a of the order of 40%; any larger variance in the perturbations 
produced negative value of a. The differences between the reconstructions 
and the true path remain within the acceptable range of errors. These graphs 
show that the filter has little sensitivity to perturbations (we did not calculate 
statistics here because the insensitivity holds for each individual run). 

We now estimate the parameter a from data. The filter needs an esti- 
mate of a to function, call this estimate o" a ssumed- If ^assumed ^ cr, the other 
assumptions used to produce the data set (e.g. independence of the displace- 
ments and of the observations) are also false, and all one has to do is detect 
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, (3) 


* (1) boat path 


* (2) its reconstruction 




— (3) its reconstruction (i) 




(4) its reconstruction (ii) . 


\k\ (1) 




N^----.». 


(4) 




L 











Figure 1: Some boat trajectories (explained in the text) 



the fallacy. We do it by picking a trajectory of a particle and computing the 
quantity 

D _ (EjjdX^ - dXi)f + (Ea (<*r j ' +1 - dYi)f 
J2 J 2 (dXi +1 - dXi) 2 + ^2 (dYi +1 - dYi) 2 

If the increments are independent then on the average D = 1; we will try 
to find the real a by finding a value of cr assume d for which this happens. We 
chose J = 40 (the early part of a trajectory is less noisy than the later parts). 

As we already know, a single run cannot provide an accurate estimate of 
a, and accuracy in the reconstruction depends on how many runs are used. 
In Table III we display some values of D averaged over 200 and over 5000 
runs as a function of the ratio of cx assumo d to the value of a used to generate 
the data. From the longer computation one can find the correct value of a 
with an error of about 3%, while with 200 runs the uncertainty is about 10%. 

Table III 

The mean of the discriminant D as a function of cr assume d/o", 30 particles 
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rr i / rr 

u assumed / u 


^finn rnrm 

JUuu 1 Lillo 


onn Time; 




J- . J_ t: 1 . W _L 


J_ . _ JL _l_ . WO 


6 

. W 


1 08 + 01 

± . WO _l_ . W -L 


114+ 07 

J. . -L t: _1_ . W 1 


7 


1 05 + 01 

J. . WtJ _l_ . W -L 


110+ 07 

J. . J. W _l_ . W 1 


.8 


1.04 ± .01 


1.14 ± .07 


.9 


1.00 ± .01 


1.01 ± .07 


1.0 


1.00 ± .01 


.96 ± .07 


1.1 


.97 ± .01 


1.01 ± .07 


1.2 


.94 ± .01 


.99 ± .07 


1.3 


.93 ± .01 


1.02 ± .07 


1.4 


.90 ± .01 


.85 ± .06 


1.5 


.89 ± .01 


.93 ± .07 


2.0 


.86 ± .01 


.78 ± .05 



7 Conclusions. 

We have exhibited a non-Bayesian filtering method, related to recent work on 
chainless sampling, designed to focus particle paths more sharply and thus 
require fewer of them, at the cost of an added complexity in the evaluation of 
each path. The main features of the algorithm are a representation of a new 
pdf by means of a set of functions of Gaussian variables and a resampling 
based on normalization factors. The construction was demonstrated on a 
standard ill-conditioned test problem. Further applications will be published 
elsewhere. 
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